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Inspiral of binary black holes occurs over a time-scale of many orbits, far longer 
than the dynamical time-scale of the individual black holes. Explicit evolutions 
Q> ' of a binary system therefore require excessively many time-steps to capture inter- 

■ esting dynamics. We present a strategy to overcome the Courant-Friedrichs-Lewy 
psj . condition in such evolutions, one relying on modern implicit-explicit ODE solvers 

, and multidomain spectral methods for elliptic equations. Our analysis considers the 

, ^ \ model problem of a forced scalar field propagating on a generic curved background. 

I Nevertheless, we encounter and address a number of issues pertinent to the binary 

' black hole problem in full general relativity. Specializing to the Schwarzschild ge- 

ometry in Kerr-Schild coordinates, we document the results of several numerical 
O . experiments testing our strategy. 
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I. INTRODUCTION 



> 

0^ ■ Numerical simulations of the inspiral and merger of binary black holes (BBH) investi- 

^ . gate Einstein's equations in the nonlinear regime where analytical progress often proves 

\ intractable. The primary goal of these simulations is the computation of gravitational wave- 

O ' forms necessary to analyze output from gravitational wave detectors like the "Laser Interfer- 



ometric Gravitational Wave Observatory" (LIGO). Breakthroughs in 2005 have yielded two 
ways to simulate BBH evolutions: the generalized harmonic system (GHS) with excision hj 
and the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) system with moving punctures [J,!^. 
^ ' Over the last few years numerical relativity has seen rapid progress along both fronts. 

■ The evolution of a binary black hole proceeds through three phases. During the inspiral 

phase, the two separate black holes orbit about each other, with the orbit gradually tighten- 
ing due to emission of angular momentum and energy via gravitational radiation. At small 
separation, the black holes encounter a dynamical instability, plunge rapidly toward each 
other and merge. This merger phase results in a single, larger, highly distorted black hole 
which subsequently relaxes to a stationary black hole during the ringdown phase. Merger 
and ringdown happen quickly, together lasting about 200M, where the black hole mass M 
sets both the spatial and temporal scales. Therefore, merger and ringdown are compara- 
tively easy to simulate at modest computational cost. In contrast, simulation of the inspiral 
phase is a daunting computational challenge. Because the orbital period increases rapidly 
with separation of the black holes, simulation of even a modest number of orbits requires 
much longer evolutions. For example, the last 10 orbits of an equal mass non-spinning bi- 
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nary black hole last about 2000M, already an order of magnitude longer than merger and 
ringdown. Beyond necessarily longer time-spans, inspiral simulations also require higher ac- 
curacy. Indeed, gravitational wave flux decreases with separation, and it must be accurately 
resolved in order to compute the correct phasing of the gravitational waves. 

To date all binary black hole simulations have employed explicit time-stepping, generally 
the method of lines with an explicit ODE scheme like the classical fourth-order Runge- 
Kutta method. Without question explicit time-stepping is appropriate for both merger and 
ringdown. However, during the inspiral phase, the relevant physical time-scale on which 
the binary separation changes is much longer than the dynamical time-scale M of each 
black hole. Nevertheless, the Courant condition associated with an explicit time-stepper 
heuristically requires that time-steps are proportional to the smallest grid spacing, and 
therefore explicit binary evolutions use time-steps that are typically of the order M/lOO 
to M/10. For instance, a recent 16 orbit simulation [3] required nearly 200,000 explicit 
time-steps. This issue becomes more pronounced when modeling black holes with unequal 
masses Mi > M2. The orbital period is proportional to the total mass M = Mi + M2, 
whereas the Courant limit dictates that the time-step is proportional to the smaller mass 
M2. The number of explicit time-steps needed to ensure numerical stability then scales like 
M/M2. Due to these reasons, only a few binary black hole simulations with mass ratios above 
4:1 have been performed jl, [6|, and these are quite short and computationally expensive. 
Courant limitations are likewise more severe for simulations of spinning black holes, which 
also require higher spatial resolution close to the black holes. 

These arguments suggest that some form of implicit time-stepping would more efficiently 
treat the inspiral phase, and this paper begins the study of alternative ways to carry out 
temporal integration of orbiting binaries in the early phase of their evolution. Our approach 
is based on modern implicit-explicit (IMEX) ODE solvers 0, B 0] and classical multido- 
main spectral methods^, and in particular those [ll| used for solving elliptic problems (it 
turns out that our implicit equations correspond to elliptic PDE). The generalized harmonic 



formulation [12|, |13|, |lj] rewrites Einstein's equations as 10 scalar wave equations for the 



components of the metric, which are coupled through nonlinear lower order terms. Because 
the principal part is just the scalar wave operator on a curved background spacetime, we 
consider here the model problem of evolving a scalar field on a single black-hole background. 
IMEX schemes like the ones pursued here are not the only possible approach to circum- 
vent the Courant-Friedrichs-Lewy condition. For instance, Hennig and Ansorg [15] explore 
spacetime spectral methods to solve scalar wave equations. 

The organization of the paper is as follows. The upcoming Section[TTlgives a brief overview 
of IMEX methods, using the specific example of Additive Runge-Kutta (ARK). It also briefiy 
collects the relevant first-order equations describing the propagation of scalar waves on a 
generic curved spacetime, and discusses boundary conditions for such equations. Section UTTl 
contains our main analytical discussion, and it focuses on the novelty of solving the implicit 
equations which arise in our time-stepping strategy. Much of this theoretical analysis is 
general, but we eventually settle on the concrete example of a scalar field propagating on 
the Schwarzschild geometry in Kerr-Schild coordinates. Section IIVI describes the results of 
several numerical experiments carried out for the Schwarzschild scenario. The conclusion 
in Sec. |V] summarizes our findings and discusses steps necessary for application of IMEX 
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methods to the ultimate target problem, binary black hole inspiral. Finally, three appendices 
collect some technical calculations omitted in Section UTTl 

II. PRELIMINARIES 
A. Implicit-explicit additive Runge-Kutta 

From the computational point of view, all IMEX methods require that we are able to 
numerically solve an implicit equation. For concreteness, we here consider ARKS (2) and 
ARK4(3), two IMEX additive Runge-Kutta schemes introduced in These schemes share 
the same algorithmic structure (only their sets of Butcher tableaux differ). ARK3(2) is 
a 4-stage third-order scheme with a second-order embedded scheme, while ARK4(3) is a 
6-stage fourth-order scheme with a third-order embedded scheme. Although we will not 
report on it here, we have also considered various versions of semi-implicit spectral deferred 
corrections (SISDC) 0, Q as an alternative to ARK. The nature of the SISDC algorithm is 
quite different, but its implementation also requires that we are able to solve (in this case 
at each substep) an implicit equation of the same form. 

We will not discuss accuracy and stability properties of ARK. Our purpose here is simply 
to describe the algorithm, highlight what is needed for implementation, and focus on the 
origin and structure of the implicit equation. When considering first-order systems for scalar 
wave propagation below, we will adopt what is essentially a reversed semidiscrete picture. 
That is to say, we consider time as discrete, but retain the spatial continuum. When adopt- 
ing that picture, we will write down a continuum implicit equation (a spatial differential 
equation) that corresponds to the implicit equation appearing in the ARK algorithm. Al- 
though this should not prove cause for confusion, we have nevertheless raised this issue now, 
since we adopt a similar notation whether or not the spatial continuum is retained. 

Mostly adopting the notation of 0], we begin with a generic initial value problem 

'^ = f{t,u) = Y,f^\t,u), u{t,) = u,, (1) 

with u a vector of unknowns. Adopting a 2-additive scheme, we have split the right-hand 
side / into explicit (nonstiff) = /'^l and implicit (stiff) = /^^^ sectors. The ARK 
schemes specify a rule for advancing the vector n" at a present time-step t"^ (perhaps to) to 
the vector n""*"^ at the next time-step t""*"^ = + At, and this rule requires the construction 
of s stage values u^^\ z = 1, 2, . . . , s, corresponding to intermediate times t^*^ = t"'+CiAt. The 
first stage is given by u^^^ = u", and the remaining stage values are determined sequentially 
by 

i 

ii» = + AtJ2 [af^f^{t^'\u^'^) + a{^f{t^^\ n(^'))] , 2 < ^ < s. (2) 
i=i 

After all the stages have been computed, the updated solution is given by the stage expansion 

s 

^„+l = + At ^ 6i [/^(t», w») + /^(t«, w»)] . (3) 

i=l 
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The parameter matrices = (a^) and A^^°^^^ = (oL), along with the coefficients b = 

and c = (q), stem from Butcher tableaux collected in ERK stands for explicit Runge- 
Kutta, and ESDIRK for explicit singly diagonally implicit Runge-Kutta. In the ESDIRK 
acronym, the explicit refers to the trivial ffist stage, and diagonally implicit to the fact that 
the sum in ([2]) stops at i rather than s. 

For the explicit sector, = for j > i. Therefore, in Eq. ([2]) the ffist term in the 
sum depends solely on already known stages, u^^\ . . . ,u^'^~^\ In contrast, for the implicit 
sector, a^j = for j > i, with af^ = 7 7^ unless i = 1 (the singly in ESDIRK indicates 
that the diagonal elements = ^33 = ■ ■ ■ = al^ all equal the same constant 7). The term 
a^if\t^'^\u^'^^) turns Eq. ([2]) into an implicit equation for u^'^\ Implementation of an ARK 
scheme therefore requires that we are able to solve (at each stage after the ffist) an implicit 
equation of form 

u-af\t,u) = B, (4) 
where a = 'yAt and B depends on the previous stage values. 



B. First-order equations for a scalar field on a curved background 

Our goal is to solve the scalar wave equation 

V^V^t/; = S (5) 

with a given source term S = S{t,x''). We consider Eq. ([5]) on a generic curved background 
with line-element given in the usual 3+1 decomposition, 

ds^ = -N^dt^ + gjk {dx^ + VHt) {dx'' + V^dt) . (6) 

Here gjk is the induced metric on t = const hypersurfaces, is the lapse function, and 
is the shift vector. These quantities are known functions of space x^ and time t, where 

lower-case Latin indices (7, fc, . . .) denote spatial components running over 1, 2, 3. 
Following Hoist et al. [l6|, we rewrite Eq. ([5]) as the following ffist-order system: 

d^ip = V^dkip - NH (7a) 
dtH = V'^dkU - Ng^'^dj^k + NKU + N^kJ'" + NS (7b) 
dt^j = V^dk^j - Ndjil + ^kdjV^ - UdjN. (7c) 

Apart from the possible inhomogeneous forcing term NS in ( ]7bll . these are Eqs. (14-16) of 
n is a new evolved variable representing the time derivative of ip, and Eq. (|7aj) is the 



definition of H. The $j = djip represent the spatial derivatives of ^p. The quantities J'^ and 
K depend only on the background spacetime, 

= -N-^g-^/^dj{Ng^/^g^'') (8) 
K = -N-'g-'/'[d,g'/'-d,ig'/'V% (9) 

where K is the trace of the extrinsic curvature tensor and g = det{gjk). These formulas for 
and K are respectively Eqs. (17) and (18) of [l6l |. 

Solutions of the ffist-order system Eqs. ([7]) are equivalent to those of Eq. ([5]) only if the 
constraint 

Ck = dkiJ - (10) 
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vanishes. The constraint provides an important hnk between Eq. (17a|) and (17c|) . as the latter 
equation is derived by taking the time derivative of the constraint: 

dtCk = dt^, = dtdk^ = dkdt^. (11) 



Thus, the right-hand side of Eq. flTcl) is the gradient of the right-hand side of Eq. flTal) . with 
dkip replaced by 

Boundary conditions relative to a boundary element with outward-pointing unit normal 
Uk are described in terms of the characteristic fields 

Z' = ij, Z| = P/$fc, U'^ = U±n'^j,, (12) 

where Pj^ = gj — n^Uj indicates projection tangential to the boundary element. Relative to 
the time axis d/dt, the coordinate speeds of these fields are respectively 

-nkV\ -nkV\ -nkV^±N. (13) 

A characteristic field requires a boundary condition whenever its characteristic speed is 
negative. The scenario we consider later, that is wave propagation on a Schwarzschild black 
hole, has two boundaries: First, an outer spherical boundary Bo where Z^, Z'j and are 
incoming and require boundary conditions. Second, an inner spherical boundary Bi which 
is inside the black hole horizon and surrounds the singularity at the center of the black hole. 
On the inner boundary all characteristic fields are outgoing (i. e. moving toward the center 
of the black hole), and boundary conditions must not be imposed on it. This pure outflow 
boundary results in several interesting features of the present work to be discussed below. 
If Bo and Bi are adapted to the background symmetry (i. e. round spheres, which they need 
not be for our numerical work), then n^d/dx^ oc d/dr on Bo and n^d/dx^ oc —d/dr on Bi. 

The boundary condition on is physical; this boundary condition and the choice of 
initial data determines which solution of the second-order wave equation ([S]) is computed. In 
this paper, we typically choose the initial data, the boundary values of f/^~, and the external 
forcing S such that the solution follows a prescribed exact solution. Boundary conditions on 
the fields Z^ and if necessary, are chosen to ensure that the constraint Ck vanishes on the 
boundary. Solutions to the first-order system ([7]) which violate the constraint = are not 
admitted either by the scalar equation ([5]) or the reduced system which arises from setting 

= dk'ip in d?]). The boundary conditions on Z^ and Zj rule out these spurious solutions 
to the extended system ([7]), provided that the initial data also satisfies the constraint. Such 



constraint preserving boundary conditions have been derived by Hoist et al. [16| and refined 
by Lindblom et al. [ij] . The implicit equations that we encounter in our use of ARK methods 
require boundary conditions that parallel those of the evolution problem. We will therefore 
consider as given boundary data for the implicit problems, as well as the boundary data 
Ck = 0, whenever necessary. 

III. IMPLICIT EQUATIONS 
A. First-order equations 



The ARK algorithm described in Sec. Ill Al is applicable only to systems of ODE, and for 
the case at hand such an ODE system arises upon spatial approximation of Eqs. (JT]) via a 
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pseudospectral collocation method (see, for example, [l6| for details). However, as mentioned 
earlier, we find it convenient to retain the spatial continuum in our discussion, and so write 
down the continuum implicit equations (PDEs) which, upon spatial approximation, yield 
the relevant algebraic implicit equations appearing in our IMEX algorithms. Equations ([7]) 
have the form of Eq. ([T]) for the evolved variables u = {ip, 11, We will consider a number 
of possibilities for splitting the right-hand side of Eqs. (JTj) into stiff (implicit) and nonstiff 
(explicit) sectors, but always treat the system's principal part (i. e. all spatial derivatives) 
implicitly. 

Each of our possible choices for the IMEX splitting is specified by writing down the field 
components of the implicit equation (jl]). Treating implicitly the first two terms from each 
right-hand side in ([7]), and possibly the forcing term NS from fl7bl) . we get case (i): 

^P-a (V^dmiJ - NU) = (14a) 
U-a {V^dmll - Ng^"'dj^m + eNS) = Bn (14b) 
$fc - « {V'^dA - NdkU) = (14c) 

with e = 1 for implicit treatment of NS, and e = otherwise. Therefore, as with the other 
cases to follow, case (i) is actually two cases. A second, and similar, set of equations stems 
from also treating implicitly all terms in the right-hand side of (!7cl) . Namely, case (ii): 

^ - a{V"'dmip - NU) = B^ (15a) 
n - a{V"'dmIl - Ng^'^djl^m + eNS) = Bn (15b) 
- a{V'^dA - NdkU + ^mdkV^ - lidkN) = B^^. (15c) 

Finally, treating all or nearly all terms implicitly, we arrive at case (iii): 

^ - a{V'^d^^l) - NU) = B^ (16a) 
n - a {V"'dmU - Ng'^^dj^m + NKIi + N^mJ"^ + eNS) = Bn (16b) 
-a{V"'d.^^k - NdkU + ^.mdkV"' - UdkN) = B^^. (16c) 

For each of our three cases, we note that the inhomogeneity B = {B^, Bn, B<s>^} corresponds 
to the term in (jl]) built with ARK stage values, cf. Eq. ([2]). While we have considered only 
three possible IMEX splittings (really six including the e = 0, 1 choice), other variations are 
of course possible. Ignoring the subcases afforded by the choice of e, case (i) corresponds 
to the minimal implicit sector for which our methods are applicable, case (iii) to the fully 
implicit scenario, and case (ii) to a scenario in the middle. Note that for cases (ii) and (iii) 
the gradient of the left-hand side of the ip equation [Eqs. fll5ap and fll6ap . respectively] gives 
the corresponding $fc equation [Eqs. fll5cp and fll6cp . respectively], up to the replacement 
dk'ip — * $fc- This mirrors the structure of the first-order PDE, cf. the remark after Eq. (fTTIl . 

While we do not solve these first-order systems numerically, we expect that the follow- 
ing theoretical considerations are relevant. We view each set [Eqs. ([HD, (fT5|) . or (fT6|) ] of 
implicit equations as a spatial boundary value problem subject to Dirichlet boundary con- 
ditions on the same characteristic fields fll2p as those described in the last paragraph of 
III B[ In other words, the choice of boundary data for these implicit solves corresponds to 
the same boundary data controlled in the evolution initial-boundary-value problem. This 
physically reasonable viewpoint is analyzed further in an appendix. On the outer bound- 
ary Bo where V'^Uk > 0, we fix Z^, Zj, and as boundary data. Typically, V'^Uk < 
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on Bi, so no boundary conditions on and Zj are imposed. If —V'^ni^ — N < on Bi, 
then we would fix f/^~ as boundary data. This inequahty would hold, for example, in our 
Schwarzschild scenario, provided Bi were chosen as a surface outside the horizon. Now, each 
set [Eqs. dHl), ( fT5l) . or (fT6l) ] involves the first-order derivatives of 5 fields, whence we expect 
that 5 boundary conditions are needed to uniquely determine a solution. Indeed, 11 and 
$fe should be determined by ffT^.c). ffTSb.c) or f|T6b.c) and specification of the following 4 
boundary conditions: Zj and on Bo, and on Bi (provided —rikV'' — N < 0). Once 
n is known, the remaining equation for ip could then be integrated subject to a remaining 
fifth boundary condition for Z^ on Bo- We analyze a simplified system which justifies this 
counting argument in Appendix |Al 

For the Schwarzschild scenario when Bi lies inside the horizon, the situation is different. 
Let us view Bi as spherically symmetric, and let us extend the normal n'^ to Bi smoothly 
into the volume such that n'^ is normal to r = const spheres, and normalized such that 
Qij'rtn^ = 1. Combination of the first-order implicit equations for 11 and yields 

U^- - a [{V'' + Nn^)dkU^- + ..■]= 5n - n'^B^^, (17) 

where we define U^^ = 11 — n^^k even away from Bi. On the horizon + Nn^ = 0, and 
thus Eq. (fTTI) determines algebraically. Integration of Eq. (fT7|) inward from the horizon 
to Bi then results in the value of on Bi. Thus on Bi is determined self-consistently 
by the equations, and we are not free to pick it. 

As with the evolution initial-boundary-value problem, in our implicit boundary value 
problems we relate some boundary data to the constraint C^. First, we identify the tangential 
components Pj^Cklso with the boundary data Zj\q^. In other words, on Bo we set Z| = 
Pjidk'ip — Ck), where PjCk is a fixed function (typically zero). Along with the boundary 
data f/^~, these tangential components then allow for recovery of 11 and ^k- Second, writing 
f im . ffTSal) . or ffT6all as 

iP - aV'^Ck = aiV'^^k - NU) + B^, (18) 

we may now view all terms on the right-hand side as a given source. Rather than fixing 
Z^ = xp as boundary data on Bo, we equivalently fix V'^Ck\j3^, since ip\t3^ can then be recovered 
from the last expression evaluated at Bo- We may then formally view our outer boundary 
conditions on Bo as controlling and C^. 

For cases (ii) and (iii) the listed implicit equations determine an implicit equation for the 
constraint. Indeed, notice that the pairs fll5ap . fll5cp and fll6al) . fll6cl) are the same. If we 
subtract, say, fll6cp from the Cartesian derivative of fll6ap . then we arrive at 

Cu - a{V'^d^mCk + CmduV^) = d^B^ - (19) 
an equation we may alternatively express in terms of the Lie derivative as 

Ck — a£vCk = dkB^ — B^^. (20) 
Contraction of f[T^ on aV^ yields 

aV% - a^V^djiV^'Ck) = aV^d^B^ - 5$ J. (21) 

In principle, Eq. fl20|) might be integrated along the shift, say inward from the outer boundary 
Bo where Dirichlet boundary conditions on Ck are set. 
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The ARK scheme should preserve the constraint; i.e. if = initially, it should remain 
zero. We investigate this point by combining ([2]) and (jl]) into a formula for the ith stage 
source, 

i-l 

B« = li" + At ^ [afjfit^'K u'^'^) + aijf{t^'\ u^'^)] , 2<t<s, (22) 
i=i 

where we recall that u^^^ = u^. For both cases (ii) and (iii), the if) and $fc components of 
f^{t,u) vanish, whereas 

/ V'dkilJ -NU \ 

fit,u)=l . , (23) 

\ - NdkU + ^mdkV"' - UdkN J 

with • indicating an expression irrelevant for the present discussion. Therefore, if the previ- 
ous stage values u^^\ j = 1, . . . , i — 1, obey the constraint (fTOjl . that is C^^ = ■ ■ ■ = C^^'^^ = 0, 
then the ith. source will satisfy dkB^ = 3^. As a result, ([19]) will be a homo geneous 
equation for the constraint C^*'' at the zth stage, with the solution C^*^ = in the interior 

because Cf^ = has been enforced on the boundary. We will draw on these observations 
below. 



B. Second-order implicit equation 

In principle, one could solve directly the first-order implicit equations given in Eqs. dHD, 
f[T^ . or flTB]) . and we have done so in spherical symmetry. However, for the more demanding 
3d cases leading toward our ultimate goal of handling binary black holes, we would like to use 
the multidomain spectral EllipticSolver which is part of the Spectral Einstein Code 
SpEC used for binary black hole evolutions j3, E, 3]- The EllipticSolver has been written 



to handle second-order elliptic equations. Moreover, preconditioning strategies for second- 
order elliptic equations are well understood relative to those for first-order equations. For 
these reasons, we have chosen not to directly solve first-order equations. Rather, we first 
solve a single second-order scalar equation for ip, one stemming from combination of the 
above equations and subject to appropriate boundary conditions discussed below. This ip 
equation is different for each of the three cases, and it is only for cases (ii) and (iii) that we 
can show, at least formally, that our solution process is consistent with solving the original 
first-order set of equations. Once a solution ip has been determined for each case, we obtain 
n algebraically using fll4al) . (11 Sal) , or (I16ap . all the same equation. Finally, we recover 
from ip via differentiation. Therefore, at each stage in our IMEX algorithms we perform 
what amounts to a naive constraint projection. This is necessary for case (i), but would 
seem not strictly necessary for cases (ii) and (iii). 

Both ARK3 and ARK4 have explicit first stages, for which no implicit solve needs to be 
done. Nevertheless, in order to achieve stability in our IMEX evolutions, we must perform 
the naive constraint projection on the first-stage fields, at least when such projection is 
carried out on the other stages. The discussion after Eq. fl2^ pertains to exact arithmetic, 
whereas round off errors in the stage expansion ([3]) will result in a u""''^ which violates the 
constraints. While this violation is negligible over a single time-step, such violations appear 
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to accumulate. Projection of the first-stage fields guarantees that dkB^ = B\,'^ throughout 
the evolution. 

Let us first derive the scalar equation for ip associated with case (i), Eqs. f|T4l) . Combi- 
nation of (fT^ .b) eliminates the term proportional to 11, 

iP - aV^dm^ + a^NV'^d.mIl-a^N'^g^""dj^m = - aNBu - ea^N'^S, (24) 

whereas from (I14cl) we obtain 

aV^^m - a^V^V'^dj<^m + a^W^^^n = aV^B^^. (25) 

The difference of the last two equations is independent of 11, 

iP - aV^idm^ + $„) - {N^g^"" - V^V"") dj<^m = B^ - aNBu - aV^'B^^ - ea^N^S. 

(26) 

Next, we use the constraint ffTOl) to replace $fc by d^ip — C^, and find 

ijj-2aV^djtlj - [N^g^^ - V^V^) djdkip = 

B^ - aNBn - aV^B^^ - ea^N'^S - aV^Cj - a^{N'^g^^ - V^V^)dfk, (27) 

which is our second-order ip equation for case (i). 

Derivation of the ip equation for cases (ii) and (iii) is more complicated, but nevertheless 
follows the same steps taken for case (i). For example, Eqs. (ITSk.b) again lead to (jUj). 
Similar to before, we eliminate the term a^NV^dmJl with the contraction of aV^ on the 
$fc equation fllScI) . However, now this third equation is more complicated and features an 
extra factor of IT. Therefore, we first use (11 Sal) to rewrite (I15cp as 

- a,^ - a(y™9^$fc - Ndkli + ^mdkV^ - akV^d^i^) = B^^ - a^B^, (28) 

where ak = dk log A^. Finally, we contract the last equation on aV'^, subtract the result from 
flM|) . and then use ffTOj) to replace all $fc terms by dkip — Ck- These steps yield the following 
equation for case (ii): 

(1 + aV''ak)tp - [2aV'' + c?[a^V'V^ - V'djV^)\dki) - a\N^g^^ - V'V^)djdkiJ 
= (1 + aV^ak)B^ - aNBu - aV^B^^ - ea^N^S 
- a^N^g^^dfk + a^V^d^iV^Ck) - aV'^Ck- (29) 

Even more involved calculations using Eqs. (fT6|l similarly yield 

[1 + a(vV - NK)]ij - [2aV^ + a^{ajV^V'' - V^djV^ - N^j'' - V''NK)]dkip 
-a^{N'^g^'' - V^V'')djdktp = [l + a{V^ak - NK)]B^ - aNBu - aV'^B^^ 

- ea^N^S + a^N^{J^Ck - g^^dfk) 
+ a^V^d,{V''Cu)-aV^Ck, (30) 

the second-order ip equation for case (iii). Notice that neither fl27|) . fl29|) . nor fl30|) features 
derivatives oi B = (i?^, i?n, -B$s,)- The absence of B derivatives indicates that we have not 
differentiated any of our original first-order equations. 
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Let us now consider the corresponding boundary conditions for the second-order equations 
( !27|) . (!29l) . (!30|) . Combination of fll4al) . the same for all cases, with the formula (fT2|) for ?7^~ 
yields 

^ - aV'^dm'^ + aN{U^- + n'^^fc) = 5^. (31) 
Using the constraint = dkip — ^k, we therefore find 

^ + a {Nn"" - 1/™) dmip = B^- aNU^' + aNn'^Ck (32) 

as our boundary condition. If Bi lies inside the horizon, we shall not impose ( !32l) on as the 
initial boundary value problem for Eqs. ((71) does not require an inner boundary condition. 
For a Schwarzschild spacetime in Kerr-Schild coordinates. Sec. IIII CI shows that an outer 
boundary condition alone is indeed sufficient. In this setting, the inner boundary condition 
for the second-order equation is replaced by a regularity condition across the horizon. 
Numerical tests in Sec. IIVI indicate that this viewpoint holds in a more general setting. 
While we have been careful to retain all constraint terms in deriving Eqs. (1271) . (!29l) . (!30|) 
and (132|) . in practice we have set = before solving these equations numerically. 

At least for cases (ii) and (iii), the following argument formally proves that — in lieu of 
directly solving the full set of first-order implicit equations (fT5l) or (fT6|) — we may instead take 
the following steps: solve (1201) . solve either (12^ or (150]) (the single second-order ifj equation), 
and then recover 11 from Eq. fllSap or fll6al) and $fc from Eq. fllUp . In combining the first- 
order equations to produce the single second-order tp equation, we have not differentiated 
the original equations, rather second derivative terms have arisen via substitution with the 
constraint dj(^k djdkip — djCk- Therefore, the ip equation is truly a second-order equation 
only if we interpret the constraint terms appearing on the right-hand side as part of the 
inhomogeneity. To achieve this interpretation, we assume that we may integrate the system 
fl20|) . inward from the outer boundary Bo, where we fix Cklso- These boundary conditions 
are consistent, as we have earlier argued that Ck\Bo P^'^t of the boundary data for the 
original first-order set of implicit equations. 



C. Schwarzschild geometry 

We now specialize the above equations to the Schwarzschild geometry written in Kerr- 
Schild coordinates. The line-element is 

= -N^f + (dr + V^dt) ' + (d^^ + ^^^2 ^^^2^ ^ (33) 
where the lapse, radial lapse, and shift are given in terms of the mass parameter M by 



Vr+2M' V r ' r + 2M ^ ' 



This is the same line-element as given in Eq. (59) of [16 . 

Consider the Cartesian coordinates stemming from the polar coordinates (r, 9, 0) via 
the standard formulas: {x,y,z) = {r sin 6 cos (f),r sin 6 sin (f),r cos 6). We introduce a radial 
vector u'^ = x^/r which is not normalized with respect to the spatial metric determined 
by (I33f34l) . The vector = L~^v^ is the outward-pointing unit normal to the spherical 
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foliation of a spacelike level-t hypersurface. With respect to the Cartesian coordinates x'', 
we may express the spatial metric and inverse metric as follows: 



,jk 



-2 



(35) 



Here 6jk = = diag(l, 1, 1) is the fiat metric, and djr = 5jkV^ . To avoid ambiguity as to 
which metric {gjk or 5jk) has been used to lower the index, we will not write z/^. With the 
above formulas for the Cartesian components of the metric and Eqs. we find that 



J' 



1 

l2 



L' 
T 



N 



K 



1 



2 T' 

+ -V'' H V 

r L 



(36) 



To reach these equations, we have used the following identities (valid in the Cartesian coordi- 
nate system): g^^"^ = L, V'' = V^u'', uWjV^ = 0, and dju'' = r^^(Sj — UjV^). Equations (!35|) 
and (l36l) hold for any choice of A^, L and . Specializing to the values given in Eq. (I3i~ 
we obtain 



J'' 



2M(r + 4M) 
r(r + 2M)2 ^ ' 



NK 



2M{r + 3M) 
r(r + 2M)2 ' 



(37) 



For the chosen Schwarzschild background and coordinates, we now show that our solution 
procedure involving the second-order ip equation ( !29l) or ( l30l) is equivalent to solving the 
original set of first-order equations ( fT5i) or ( |T6l) . In establishing this claim we must show 
that (fT^ can be integrated inward from the outer boundary ;Bo at r = rmax, and that a 
reg'«/ar solution to the ip equation is determined by the outer boundary condition alone. We 
consider the integration of ( fT9i) in Appendix [Cl and turn to the latter issue now. Whereas 
integration of (fT9l) is only relevant for cases (ii) and (iii), the issue of a regular solution to 
the ip equation also pertains to (l27j) . and so we include this equation in our analysis. Each 
of the second-order scalar equations [Eq. fl27j) . (129!) . or (!30|) ] takes the following form: 



7^(r)^ + aS{r)u''dk'^ + {N'^g^'' - VW^) djdkip = Q. 



(38) 



Here we view constraint terms appearing in if present, as predetermined via integration 
of ([in]). We continue by calculating 



g^^djdkilJ = L-^d'^ip + 2r-^drip + r-^As2ip, 



(39) 



where A52 is the Laplace operator associated with S"^, the unit-radius round sphere. Next, 
we set = 4'em{r)Yim{(^, (p) in dSSD, thereby obtaining the equation for a generic spherical- 
harmonic mode, 



7^(r) 



+ 1) 



aS(r) + 



1? - 



This equation has the form 

Q{r)w + aV{r)w' + a^{r - 2M)w" = h{r) 

with 

a'^N^l{i + 1) 



Q{r) = (r + 2M) 



n{r) - 



-p(r) = {r + 2M) 



S{r) + 



2aiV2 



(40) 
(41) 

(42) 
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FIG. 1: Field configurations at initial and final times. 



Note that the coefficient of the second-order term w" passes through zero at the black hole 
horizon, r = 2M. We study Eq. ( HTi) in Appendix [HI where we show that r = 2M is a 
regular singular point. Furthermore, in the appendix we compute the indicial exponents 
associated with the singular point, and argue that, despite the second-order character of 
(HTj) . an outer boundary condition alone determines a unique solution which is regular up 
to and even across the horizon. 

Thus the following picture emerges: If the radius r^i^ of Bi satisfies rmin > 2M, then the 
scalar wave equation requires a boundary condition for on both Bi and Bo', the corre- 
sponding second-order implicit equation (155]) is everywhere regular and requires boundary 
conditions on both boundaries as well. For rmin < 2M, the inner boundary is an outfiow 
boundary for the scalar wave equation, and a boundary condition on is necessary only on 
Bo] in this case, a unique solution to Eq. (155]) is determined by an outer boundary condition 
alone, along with the assumption that the solution is regular across the horizon. 



IV. NUMERICAL TESTS 
A. Comparison with explicit time-stepping 

In this first subsection we demonstrate that our numerical IMEX algorithm can solve a 
standard initial value problem. To do so, we consider the evolution of pulse initial data, a 
problem for which an explicit algorithm would be better suited. Here we are evolving the 
wave equation V^V^V = without a source term, that is for S = 0. 

The following experiment has been carried out both with a one-dimensional radial code (in 
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FIG. 2: Error of implicit evolutions relative to the explicit reference solution. The dotted line is a 
least-squares fit of the last five Id data points, although shifted to make for a better figure. 

Matlab) and with a three-dimensional code (in SpEC). Whereas the 3d code uses the vari- 
ables {V^,n, ^y, ^z}, the Id code uses {'ipyH, $ = drip}- On the radial domain [1.9, 11.9], 
initial data 

Id: V = 0, n = exp [-(r-5)2] , $ = 0, 

3d: 7/^ = 0, n = exp[-(r-5)2]Re[Fn(^,0)], = 0, (43) 

is evolved to time tgnai = 15, with the background geometry taken as M = 1 Schwarzschild 
in Kerr-Schild coordinates. Initial and final radial mode profiles are depicted in Fig. [H The 
first step of the experiment is to generate a reference solution, using an explicit Runge-Kutta 
(ERK) time-stepper, either (for Id) the classical fourth-order scheme or (for 3d) the fifth- 
order Cash-Karp scheme [l^. In both cases we choose a fixed time-step At ~ 0.00366, and 
so are not using the potential adaptivity of the Cash-Karp scheme. For both the Id and 3d 
experiments, we place no boundary condition at the inner radius r = 1.9, and a Sommerfeld 
boundary condition = [cf. Eq. (IT^ . either 11 — $ = in Id or 11 — = 

in 3d] at the outer boundary r = 11.9. We further enforce constraint-preserving outer 
boundary conditions which are analogous to the boundary conditions applied to the black 
hole evolutions in ll4l. For the scalar characteristic field and the 3d code, we use 



dtZ^ = -NU + l/'^^fc, (44) 
cf. Eq. (40) Ref. 0. For the Id code we similarly use dfZ^ = — A^II + V^^ as the outer 



boundary condition. For Zf, we employ the analogue of Eq. (65) in Ref. |1^ 

dtZf = DtZ^ - nkV'n"'{dm^, - 9,$^). (45) 
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FIG. 3: Performance of ARKS for case (i) IMEX splitting. The dotted line corresponds to exact 
third-order convergence. 

Here DtZ^ denotes the P/ projection of the right-hand side of Eq. (!7c|) . We use one spherical 
shell, with 61 radial collocation points. We fix the angular resolution for the 3d evolution 
with £max = 5 as the top spherical-harmonic index. The explicit evolution uses the same 
angular filtering as [l6|. When the right-hand sides of Eqs. ([7]) are computed, they are 
transformed to scalar spherical harmonics (for the ip and 11 components) or vector spherical 
harmonics (for the $fc component), and the top two modes are truncated. 

We next carry out the same evolution via IMEX evolutions taking case (iii) (the e value 
is irrelevant because S" = 0), and check the results against the reference solution. During 
an evolution, the requisite implicit solves have been carried out with the EllipticSolver 
in SpEC, as described in [ll[. For this simple class of problems, we have chosen finite- 
difference preconditioning. Precisely, we have used an exact LU decomposition of a finite- 
difference approximation Apu the operator associated with the second-order ip equation. 
The EllipticSolver interfaces with petsc's iterative solvers, and for the iterative linear 
solves we have used GMRES, choosing all error tolerances close to machine precision. As 
mentioned, for these experiments rmin = 1.9M < 2M, so the inner boundary lies inside of 
the horizon. Therefore, throughout our evolutions we solve the second-order ip equation with 
no boundary condition at the innermost collocation points. Rather in the relevant matrix- 
vector multiply needed for the iterative solver, only the PDE is enforced at these points. 
The elliptic equation for ifj is solved for Y^m modes with £ < ^max — 2 to avoid technical issues 
due to representation of Cartesian tensor components with scalar spherical harmonics. This 
acts as an angular filter for ip, obviating the need for further angular filtering as described 
above for the explicit evolution. 

For a number of temporal resolutions and for the ARK4 method, we show the results in 
Fig. [2J The errors plotted in the first, second, and third quadrants correspond to the Id 
code, and these errors have been computed after interpolation onto a finer uniform grid. 
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FIG. 4: Performance of ARKS for case (ii) IMEX splitting. 



The fourth quadrant plot collects results from both the Id and 3d experiments. The black 
circles correspond to the Loo errors from the Id code shown in the other plots (and are taken 
over all fields). The red squares are errors from the 3d experiment (and, again, are taken 
over all fields). Note that these errors have been computed in 3d. In both cases we see clean 
fourth-order convergence. 



B. Model problem on a black hole 

We now consider a problem for which temporal variations occur on time-scales much 
longer than the Courant limit for an explicit time-stepper. This model mimics the binary 
black hole configuration, providing a testing ground for our IMEX methods. Our model 
problem on a single M = 1 black hole is set up as follows. For the solution we adopt the 
Ansatz 

Mt,x,y,z) = cos(c^t)/(r)Re[F2i(^,0)], (46) 
where /(r) is a radial profile. This profile is chosen as a polynomial of degree 2q, 

/(r) = ( — - — ) Ur - ri){r2 - r)Y, r2 > ri, (47) 

truncated so as to vanish whenever r lies outside the interval [ri,r2]. We typically choose 
g = 5, ri = —1, and r2 = 11.9. The computational domain covers radii r G [rmin,rinax] = 
[1.9, 11.9] C [ri,r2]. Note that r^^^^ = 1.9 is somewhat inside the black hole horizon, r = 2 
(recall that M = 1), and we therefore never apply a boundary condition at rmin. Never- 
theless, we have chosen the support of / such that / is non-zero at the inner edge rmin of 
the computational domain. However, / does vanish at the outer boundary rmax, a necessary 
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FIG. 5: Performance of ARKS for case (iii) IMEX splitting. 



requirement for avoiding boundary-driven temporal order reduction |20|]. We substitute the 
chosen ipQ into Eq. ([5]) to compute the source S. We furthermore initiahze the initial con- 
ditions for ip, n and with the Ansatz ipQ. The uo value determines the time-scale of the 
temporal variations, and we present results for three values, uj = 1,0.1,0.01. For u = 1, 
temporal and spatial scales are comparable, whereas for u = 0.01, temporal variations are 
vastly slower, so that explicit time-steppers will be limited by the Courant condition. 

Our radial expansions use Chebyshev polynomials Tfc(X) as basis functions, where we 
map X e [-1, 1] to r G [rmm, ^max] via 

r(X) = Ae^^ + C, (48) 

with C = —2 and parameters A and B chosen such that r(— 1) = r^in, t(+1) = Tmax- 
The mapping fHSl) serves two purposes: First, it increases resolution close to the black 
hole, resulting in a somewhat faster convergence rate for the spectral representation of the 
Schwarzschild background, Eq. fl551) . Second, through this mapping the expansion of the 
radial profile /(r) in Chebyshev polynomials acquires non-zero power in all radial modes. 
In contrast, the linear map r{X) = X would result in only the first 2q + 1 Chebyshev 
polynomials being excited. 

With either ARKS or ARK4 and one of the considered IMEX splittings, our experiment 
is to evolve the initial data specifying this solution, assuming that the field equations include 
the exact forcing function S{t,x,y, z). We evolve for about 10 oscillation periods, to final 
time Tfinai = 65.536/ciJ. We have chosen = 25 radial collocation points, and the angular 
grid is determined by top azimuthal index ^max = 5. 

To examine the influence of the IMEX splitting on both long and short evolutions, we 
perform numerical runs with ARKS for each of the aforementioned cases (i), (ii), and (iii). 
The results are collected in Figs. [3|, HI and [51 Some of the errors in Figs. [3] and H] correspond 
to blowup and fall outside of the plot range. Comparing these plots, we notice that for 
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Error in scalar field \\f at final time 65.536/(0 for case (iii) 
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FIG. 6: Performance of ARK4 for case (iii) IMEX splitting. Tlie dotted line corresponds to exact 
fourth-order convergence. 



uj = 1 short-time runs the accuracy is insensitive to the choice of splitting. However, for 
the small-cij, longer-time runs, the fully implicit case (iii) is advantageous in the following 
sense. As u is reduced by a factor of 10 (from 1 to 0.1, and then to 0.01), this splitting 
allows for a corresponding increase of the time-step At by the same factor of 10 without loss 
of accuracy. Fig. [6] shows results for the same case (iii) experiment, but with ARK4 rather 
than ARK3 used for time integration. The fully implicit scenario, that is case (iii) and e = 1, 
corresponds to evolving solely with the L-stable ESDIRK component of the ARK algorithm 

0. 

When this same oscillating multipole problem is evolved by the explicit fifth-order Cash- 
Karp scheme 3] the Courant limit is about A^cfl — 0.235 independent of u. For all 
splittings and for all u, the IMEX code allows for time-steps one to two orders of magnitude 
larger than the explicit code. For slow temporal variations, u = 0.01, the splitting (iii) with 
e = 1 allows for a time-step At = lOOOAtcFL (i-e. about 3 time-steps per oscillation period) 
while maintaining an accuracy of about 10^'^. 



C. OfT-center model problem on a black hole 

We now consider the following off-center Ansatz: ipo{t, x, y, z) = cos{ujt)f{\r — ro|), where 
uj = 0.01 and relative to the black hole center ro = (0.5,-0.2,0.3). The radial profile 
is determined as above with ri = — 1 and r2 = 14. The numerical domain is comprised 
of three nested spherical shells, each with center Cq = (—0.08,0.05,-0.06). Therefore, 
each shell is neither a level-r surface nor a metric sphere with respect to the background 
Kerr-Schild geometry. Relative to their common center, the shells are determined by the 
radial bounds 1.8, 5.13, 8.47, and 11.8, with coordinate radial separations computed using 
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Errors up to time 6553.6 for w = 0.01 and At = 12.4 
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FIG. 7: Three-dimensional off-center experiment with ARK4. 



the background Cartesian coordinates described before Eq. (1351) . Each shell has Nr = 
15 Chebyshev-Lobatto collocation points, with a top azimuthal index i^^x = 9 fixing the 
angular grid. The numerical code expands variables in spherical harmonics centered on Cq. 
Because Cq 7^ Tq all spherical-harmonic modes are excited in this experiment. 

For this type of exact solution (which involves exact control of a nonzero U^^ as an 
inhomogeneous boundary condition), we expect temporal order-reduction, a well-known 
pitfall of exact time-dependent boundary conditions [20]. Therefore, our purpose here is 
not to consider temporal convergence, rather to demonstrate robustness of our evolution 
procedure in a setting which mixes several issues at once: an asymmetric solution, absence 
of an inner boundary condition, and multiple domains. While the inner boundary lies within 
the horizon, the coordinate characteristic speeds vary spatially across it. Fig. [7] depicts long- 
time error histories for all fields using the fully implicit ARK4 time-stepper, that is case (iii) 
with e = 1. The plot clearly shows the fields' response to the external forcing, with the errors 
continuing to oscillate. At least for linear problems we consider here, we believe that our 
implicit evolutions are robustly stable, even in the absence of an inner boundary condition. 



D. Model problem with perturbed initial data 



The numerical experiments described in Subsections IIVBI and IIV CI deal with scenarios 
in which our IMEX integration is simply driven by an external forcing, and as a result no 
secular errors accumulate. That we are therefore able to achieve reasonable accuracy for 
large time-steps is perhaps not surprising. In this section, we provide one further test which 
combines transient, rapidly changing behavior at early times with a slowly varying solution 
at late times. This test mimics start-up effects encountered in binary blackhole simulations, 
which typically exhibit rapid transient behavior at early times when the black holes settle 
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FIG. 8: Errors for the modified model experiment. Note the dotted curve giving the 
deviation of the explicit-reference ip from the Ansatz ipQ. The large deviation at early times is 
present since the Ansatz is not the exact solution. 

down from imperfect initial data into their quasi- stationary configuration. 
We still solve the scalar wave equation with the same source term as before, 

V^V^^ = S, S = V^V^V-o, (49) 

with tpoltjx'^) given in Eq. fH^ . and here with u = 0.01. However, we now choose initial 
conditions for generating ifjit^x^) which are inconsistent with those for generating ip^it^x^). 
Specifically, we choose 

V^(0,x'=) =7/'o(0,x'^) (50a) 
n(0, x^) = no(0, x'') + G{x'') (50b) 
$(0,x^) = <l>o(0,x'=), (50c) 

where G{ru^) = exp[— (r — 5)^]Re[Y'ii(6', 0)] is the angularly modulated Gaussian wave packet 
used in Sec. IIV A| see Eq. ( H3|) . Because of the presence of G{x''), the solution to this 
evolution problem is not simply ip{t,x^) = ipo{t,x^), but rather there will be an initial 
deviation. For long evolutions, the effect of the Gaussian perturbation dies away (due to 
our dissipative radiation boundary conditions), and ipit^x^) ~ ipo(t,x'') for large t. For this 
experiment tfinai = 6553.6. 

We again work with a single, centered, spherical-shell domain and Nr = 61, ^max = 5. 
Rather than the mapping fHHj) . now we choose the identity r{X) = X. This is necessary to 
fully resolve the Gaussian at early times, or else we would require even more radial points. 
We will again generate a reference numerical solution using an explicit time-stepper, in this 
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FIG. 9: Time-step sizes for explicit-reference and implicit evolutions. 



case Dormand Prince 5 (DPS) [2l|, against which we will compare an IMEX numerical 
solution obtained with ARK4, choosing case (iii) and e = 1 so that the evolution is fully 
implicit . 

A key difference between this experiment, and the ones considered in previous subsections, 
is that we now exploit adaptive time-stepping with dense output. Both DPS and ARK4 allow 
for error control and dense output. The adaptive time-stepping, based on a proportional- 



integral controller described in [21[ , allows the IMEX method to use small time-steps during 
the initial transients, and large time-steps once the transients have died away. Dense out- 
put allows us to conveniently keep track of the error history between the explicit-reference 
and implicit numerical solutions. Throughout the course of both the explicit-reference and 
implicit evolutions, we output the solution component ip ^ill times divisible by IS. For 
both the explicit-reference and implicit evolutions we choose an initial step size At = 0.04, 
with the absolute error tolerance 10~^. 

Figure [8] depicts the history of the Loo difference between the explicit-reference and im- 
plicit ip, showing that the implicit ip maintains uniform accuracy throughout the evolution. 
Also depicted in Fig. [S]is the deviation of the explicit-reference ip relative to ipQ{t,x^). Al- 
though ipoitjx'') is not the exact solution, the figure shows that it effectively is for times 
later than t = SOO. Figure M depicts the time-step sizes taken throughout the explicit- 
reference and implicit evolutions. As seen in Fig. [HI the DPS evolution essentially runs at 
a fixed time-step At ^ O.OS near the Courant limit, and in fact this evolution took over 
1.3 X 10^ time-steps. By contrast, only 47S time-steps were taken during the implicit evo- 
lution. For the implicit evolution, the time-step size starts off small and remains near 0.1 
while the Gaussian pulse propogates off the domain. However, at later times the step size 
dynamically relaxes to a time-step of order ~ 20. 
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V. CONCLUSION 

As noted in the introduction, implicit or IMEX time-stepping is bound to offer a more 
efficient means of carrying out BBH evolutions, especially evolutions involving black holes 
with markedly unequal masses. In the context of scalar waves on a single black hole, this 
paper has analyzed several issues pertinent to the eventual use of IMEX methods in ac- 
tual BBH evolutions based on the generalized harmonic system. These include the role of 
constraints, the need for second-order implicit solves, and the nature of the IMEX splitting. 

Specifically, we have investigated the role of a pure outflow boundary within the black 
hole horizon. Consistent with the physics, the initial boundary value problem associated 
with the hyperbolic system of PDEs does not require a boundary condition on the inner 
boundary. Naively, one would expect that a second-order implicit equation, as used in our 
work, would require both outer and inner boundary conditions, in disagreement with both 
the underlying physics and hyperbolic PDE. However, the second-order equation is singular 
at the horizon, and by requiring that the solution is regular across the horizon, we have 
found that the outer boundary condition alone yields uniqueness. 

We have examined the impact of different IMEX splittings on the evolution of a wave 
equation, finding that the performance of our IMEX schemes depends sensitively on the 
precise splitting choice [cf. Figs. [3] to [5]. Only the fully implicit choice [case (iii), e = 1] 
allows for time-steps proportional to the temporal time-scale, i.e. At oc 1/uj, while retaining 
accuracy independent of uj. We explain this result as follows. For small uj, the right-hand 
sides of the evolution equations ([7]) are 0{uj) by construction. However, individual terms in 
the expressions are 0{1), because tp, H and $jt are all 0{1). Only the sum of all terms is 
0{uj). Consequently, for cases (i) and (ii), and for case (iii) with e = 0, both the implicit 
and explicit sectors are large, while their sum is small. That is, and are 0{1), 
while + = 0{lj). Therefore, in the ARK scheme, both the implicit and the explicit 
sectors involve large drivings, which seems to degrade performance. In the general case, we 
conjecture that the IMEX splitting should ideally ensure that both and remain small. 

These observations suggest that a fully implicit treatment of the GHS equations will afford 
accurate evolutions with very large time-steps. However, a corresponding gain in efficiency 
may well be offset by the complexity of solving complicated nonlinear implicit equations. A 
more workable approach might be to linearize the GHS equations about the solution at the 
current time-step, in order to treat terms with constant or linear time- dependence implicitly, 
and to treat terms with quadratic (or higher) time-dependence explicitly. In this case 
would be (9(At^), and perhaps sufficiently small for rather large At. 

Another possibility for the IMEX splitting is particularly promising. Namely, splitting 
by location (or subdomain), as described in Ref. [lO(] for fluid flow past a nozzle, a problem 
for which explicit numerical evolutions are hampered by boundary induced stiffness. To 
understand the idea behind this possibility, consider the type of multidomain BBH evolu- 
tions now being carried out by the Caltech- Cornell collaboration. Such evolutions involve a 
computational domain which is split into about 60 subdomains (typically spherical shells, 
cylindrical shells, and full cylinders with axes). Among these are several concentric spherical 
coordinate shells which enclose each of the individual black holes. For either black hole, the 
innermost of these shells contains a topologically spherical apparent horizon. As these shells 
are closest to the black holes where field gradients and nonlinearities are the strongest, they 
require high resolution. Whence these shells determine the Courant limit for current BBH 
evolutions based on the generalized harmonic system with spectral methods. 
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In those shells nearest the black holes, the splitting we plan to investigate would put the 
local representation of the GHS system into the implicit sector, while the equations on all 
other subdomains would be retained in the explicit sector. The resulting evolution scheme 
would still be subject to a (milder) Courant limit arising from the grid spacing in those 
subdomains treated explicitly. However, this Courant limit would be independent of the 
resolution close to the black holes, promising efficiency gains as the mass ratio increases. 
Implicit equations would need to be solved only in a set of concentric spherical shells, 
rather than in a complicated overlapping domain decomposition, simplifying preconditioning 
and improving the efficiency of the elliptic solver. Another reason further motivates our 
interest in an IMEX splitting by location. For BBH evolutions based on the GHS system, 
implementation of outer boundary conditions (relevant only for the outermost spherical 
shell enclosing the collection of all inner subdomains) involves second derivatives of the 
physical fields [3] • The IMEX splitting we envision would treat the outermost spherical shell 
explicitly, thereby leaving in place the current implementation of outer boundary conditions. 

Finally, we point out that for black hole binaries our IMEX time-stepping strategy will 
only apply in co-rotating coordinates. Only in such coordinates does the binary configuration 
appear approximately time-independent, as the black holes remain at the same location 
in the computational grid. Moreover, the pattern of the emitted gravitational radiation 
will be almost time- independent, varying only on the inspiral time-scale together with the 
orbital frequency and the gravitational wavelength. This restriction is not onerous for our 
approach, as the SpEC code already uses co-rotating coordinates within the dual-coordinate 
frame approach developed in Ref. [isj^. 
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APPENDIX A: BOUNDARY CONDITIONS FOR FIRST-ORDER IMPLICIT 

EQUATIONS 

This appendix considers a first-order system similar to all three of our first-order implicit 
systems [Eqs. f lT^ . (fT^ . and (fTBjl ]. showing that the system requires 5 boundary conditions. 



^ In this approach inertial-frame components of tensors are evolved, and these components vary on the 
orbital time-scale. Accuracy considerations will then limit the achievable time-step to the order of the 
orbital time-scale, rather than the longer inspiral time-scale. The orbital time-scale is still a tremendous 
improvement over current time-step limitations. 
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Any of our original systems [Eqs. (fT^ . (fT5|) . or (fT6l) ] could be analyzed in a similar fashion, 
although doing so would require a mode decomposition based on vector spherical harmonics. 
Here we use simple Fourier series. Consider the system 

^ - aiV^d^tlj -U) = B^ (Ala) 

n - aiV^d^U - = Bn (Alb) 

- a{V^d,^k - dkU) = B^^. (Ale) 

where the constant shift obeys < V^^ < 1. Take the rectangular computational domain 
to be periodic in the y and z directions, and lying between x = and x = 1. Fourier 
transformation in y and z yields the transformed system 

i> - aiV^d.,ij - n) 

n - a{V^d,Il - d^, - ik2^2 - ^k3^3) 

^^ - «(\/"9ji - d^fl) 

^2 - «(^"5,,<l>2 - tk2fi) 

$3 - aiV^dAs - tksfl) 

where k2 and ks are the integers dual to y and z. All of the hatted variables should also 
carry these integer indices, e. g. 11 = n(fc2, k^), but we suppress this dependence throughout. 
We replace equations flA2bp and (lA2cp with the lightlike combinations 

U+ -a[{V^ - l)d,Ut - tk2^2 - ik3^3\ =Bn + B^, 
U- ~a[{V^ + l)d,U- - tk2^2 - ^k3^ = Bu - B^„ 
thereby arriving at the following inhomogeneous linear system: 
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With Q = ^ 
matrix A are 
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The (/c2, ^3) = (0, 0) limits of Eqs. (lAGD and (1A7P are easily computed with the result Q ~ a, 
|k| — > 0"*". The results agree with those obtained by first setting (^2,^3) = (0,0) in (lASp . 
and then performing the eigen-decomposition. Notice that the |k| = eigenvalues, which 
happen to be the diagonal entries of the coefficient matrix A in (lASp . are such that {a\g)~^ 
for g = 1, . . . , 5 are the characteristic speeds of the corresponding hyperbolic system. 
The eigenvectors ( ]A7[1 are not mutually orthogonal; however. 



4«h -«2|i^|2^^x)21g 
det [Vi, V2, Vs, V4, UgJ = ^l_(^x)2j 

and the eigensolutions 

yg{x) = e^^^Vg, g = l,2,...,5 (A9) 
form a fundamental set of solutions. Defining 

^{x) = [yi{x),y2{x),y3{x),yi{x),y5{x)] (AlO) 

and viewing the system (1A5I) as 

-^y(x) = Ay{x) + g{x), (All) 
ax 

we can now write down the general solution: 

y{x) = *(x)c + *(x) /" ^-\0s{0d^, (A12) 

where Xq is any point on the interval (0,1). The five components Cg of c correspond to 
five boundary conditions. The following recipe for fixing these components agrees with the 
convention for control of incoming fields in the corresponding evolution initial-boundary- 
value problem. The exponentials e^'^^ for q = 1,2,3,4 all blow up as x ^ 00, whereas 
gAsx (^igcays in the same limit. We want to fix the eigensolutions yq{x) for q = 1,2,3,4 
(associated with blowing-up exponentials) at x = 1, and the eigensolution 1/5 (x) (associated 
with the sole decaying exponential) at x = 0. We assume that {k2, k^) is small, so that the 
eigenvectors Vi,i^2; 1^3, and are combinations of the fields ip, U~ , $2? and $3. Therefore, 
we fix these fields at x = 1. Also for small (^2,^3), ^5 is approximately proportional to 
the fields (already fixed at x = 1) and which would be the U~ field relative to the 
outward-pointing unit normal —d/dx at x = 0. Finally then, we fix f/+ at x = 0. 



APPENDIX B: SINGULAR BOUNDARY VALUE PROBLEM 

In this appendix we consider the general solution to the second-order equation (14T|) for the 
case of the Schwarzschild geometry with respect to Kerr-Schild coordinates. Provided that 
the radial location r = rmin of inner boundary satisfies rmin < 2M, we show that a regular 
(that is, nonsingular) solution to the equation is uniquely determined by one free constant. 
We conclude that a single outer boundary conditions suffices to determine a regular solution 
to the equation. 

Our analysis assumes that both Q(r) and V{r) in Eq. fj4T]) are smooth on the radial 
domain, which is easily checked for all cases. We further note that a~^V{2M) > 1, where 
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V{2M) = 4MS{2M) + 2a. Let us verify that this last inequahty holds for the considered 
cases. For fl271) and case (i) we have 



S{r) = 2V, (Bl) 

and so a'^V{2M) = 2 + AMa'^ > 1 by ([MD- For ([29D and case (ii), we find 

S{r) = 2V' + a [{VfN-^N' - ViV)'] . (B2) 

Calculations with ([MD then show that S{2M) = 1 + |a/(4M), whence a-^V{2M) = f + 
4Ma~^ > 1. Finally, we consider ( |30l) and case (iii), with our earlier calculations giving 

S{r) = 2V' + a [{V'fN-^N' - V{V'')' - N^{J^dkr) - VNK] . (B3) 

With the formulas listed in ( !34l) and (!37l) . we find S{2M) = l — a/{4M), implying as claimed 
that a~^V{2M) = 1 + AMa^^ > 1. Our argument is completed with the following: 
Lemma: Consider the ode 

Q{r)w + aV{r)w' + o?{r - 2M)w" = h{r), (B4) 

here taken on the r-interval (2M, rmax)- Assume that Q(r), V{r), and h{r) are smooth on 
an open interval larger than (2M, rmax)- Moreover, assume that a~^V{2M) > 1, also with 
a > 0. Express the general solution as 

w{r) = CiWi{r) + C2W2{r) + wp{r), (B5) 

where wi{r) and W2{r) are solutions to the homogeneous equation (that is, for h{r) = 0), 
and Wp{r) is a particular solution. Then we may arrange for Wi{r) and Wp{r) to be regular 
as r — > 2M~^ , with W2{r) singular and obeying 

W2{r) ~ (r - 2M)^-^/("'=\ (B6) 

again as r — 2M~^. Here k = 1/V(2M), and 1 — < by assumption. The second 

solution must therefore exhibit a blowing-up (likely also branch) singularity at r = 2M. 

We begin the proof of the lemma by examining the homogeneous equation. Taken in 
standard form, that equation is 

w" + P{r)w' + Q{r)w = 0, (B7) 

where 

P(r) = Q(r) = m 

^' ar-2M' a^r-2M- ^ ^ 

Seeking solutions of Frobenius type, we then consider the indicial equation 

A(A- 1) + A/(a/t) = 0. (B9) 

Whence the indicial exponents are Ai = 0, A2 = 1 — !/(««;), and we may therefore choose 
solutions to the homogeneous problem obeying 

wi{r) ~ 1, W2{r) ~ (r - 2M)i-^/(°'^), (BIO) 
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as r — i> 2M^ . At r = 2M, the first solution is analytic, while the second exhibits blow-up, 
and likely branch behavior depending on the value of an. 

To complete the proof, we follow the method of undetermined coefficients in order to 
construct a particular solution with the desired regularity. For r > 2M, an integrating 
factor for (1B4I) is 

1 r V{0-V{2M) 



(r - 2M)-i+^/("") 



exp 



-de 



r -2M)-^+^/("''VW, (Bll) 



a J^M e - 2M 

where /x(2M) = 1. Using the integrating factor, we cast (lB4p into the following form: 

[a2(r-2M)^/("'^V(^)H'+(^-2^)"^^^^^"''V(^)Q(^)w = (r-2M)~^+^/(°^VW^(^)- (B12) 
It then follows on general grounds that 

A 

W[wi,W2]{r) = wi{r)w2{r) - W2{r)w[{r) = — — -(r - 2M)-^/("''), (B13) 

where the constant A = a'^[l — With this result for the Wronskian H^[wi,W2](r) in 

hand, we look for a solution 

Wp{r) = u{r)wi{r) + v{r)w2{r), (B14) 
subject to the variation-of-parameters Ansatz 

u'{r)wi{r) + v'{r)w2{r) = 0. (B15) 
The needed expressions for u{r) and v{r) are as follows: 

uir) = - A-'w2{0{^ - 2M)-i+i/("'^)/i(0/^(0de, (B16) 
v{r) = r A-'wriOi^ ' 2M)-'+'/^^^^ fi{0h{0dt (B17) 

J2M 

More compactly, we may write 

wpir) = f G(r,0/x(0(e - 2M)-i+i/("«)/.(0de, (B18) 

J2M 

in terms of the Green's function 

r(r n - / A''Mr)yo2{0 for 2M < r < e < 6 

'^^ \ A-iwi(0w2(r) for 2M < e < r < 6. ^ ' 

Finally, to verify that, as constructed, wp{r) remains regular as r — >■ 2M~^, we establish 
in the same limit that 

u{r) ~ Ki, v{r) ~ K2ir - 2Mf'^''^\ (B20) 

for constants Ki = u{2M) and K2 = anh{2M). The first asymptotic statement follows 
from the observation that the integrand in (]B16P is integrable at r = 2M. To get the result 
for f (r), we use 

v\r) = A-^Wi{r){r - 2M)-^+i/("''VW^(^) 

= A-^wi{2M)fi{2M)h{2M){r - 2M)-^+^/(°^) + 0((r - 2M)^/(°'^)), (B21) 

along with wi{2M) = 1 = /i(2M). Taken all together, we have shown that 

wp{r) Ki + K2ir -2M), (B22) 

as r — i> 2M~^ . Whence the lemma has been proved. □ 
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APPENDIX C: IMPLICIT CONSTRAINT EQUATIONS 

For the Schwarzschild example with hne-element fl33l) . this appendix further examines 
Eq. ( JTOl) . To obtain an orthonormal spatial triad, we complete the radial vector = L~^v^ 
defined just before Eq. ( l35l) with the standard angular directions 

eg'"' = (cos 6* cos 0, cos 6* sin 0, — sin 6*), e^*"' = (— sin 0, cos0, 0). (CI) 

In terms of the triad, we have 

Ck = UkL ^Cy + eokCe + e^kC^^ (C2) 

where = v^Ck^ Cq = ee'^Ck, and C<^ = e^'^Ck- Contraction of (fT9l) on z/'^ yields the equation 

C,-a{VX,y = B'^-u'B^^, (C3) 

with the prime denoting radial differentiation. Likewise, contraction of f|T9|) on ee^ yields 

Cg - ai^^C'e + r-'CgV') = eg'{dkB^ - 5$,), (C4) 

where in reaching this equation we have used = V^u^ and egWjU^ = r~^eg^. Similar 
manipulations establish that 

C^-a (VT; + r-^^V) = e^,' {dkB^ - 5$ J . (C5) 

Since V^'' > over the whole spherical shell, we may radially integrate (lC3|C4fC5p inward 
from the outer boundary Bo, provided C,^|b„, CqIb^, and C^\b^ are specified. We can then 
recover the Cartesian components Ck with flC2p . In principle, these components could be 
incorporated into the Q source in Eq. fl38|) . 
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